%--------------------------------------------------------------------------
% Script to generate IRFs of density differentials 
%--------------------------------------------------------------------------

clear; 
clc;
close all

set(0,'defaultTextInterpreter','latex');

% set specs 
xmin = 0;
xmax = 2.6;
xn = 100;
xgrid = linspace(xmin, xmax, xn);

nDataSel  = '1';

%--------------------------------------------------------------------------
% Load population IRFs
%--------------------------------------------------------------------------

popDir    = [pwd, '/Data/Para', nDataSel ,'/'];
pop_Dens_IRF = csvread( [popDir, 'IRF_aK_employed.csv'], 0, 0);

%%
%--------------------------------------------------------------------------
% Load IRFs
%--------------------------------------------------------------------------

nfVARSpec = '3';
nModSpec  = '3';
nMCMCSpec = '1';
nMod      = 'VAR';

sName    = ['fVAR', nfVARSpec, '_', nMod, nModSpec, '_MCMC', nMCMCSpec];
irfDir = [pwd, '/', 'Results/Para', nDataSel ,'/', sName, '/'];
PhatDens_IRF = csvread( [irfDir, sName, '_IRF_PhatDens_Aggsh1_pmean.csv'], 1, 0); % fVAR IRF
[H, ~] = size(PhatDens_IRF);
H = H-1;
Ndraws = 200;

PhatDens_IRF_uncertainty = zeros(H+1, length(xgrid), Ndraws);

for pp = 1:Ndraws
    PhatDens_IRF_uncertainty(:,:,pp)=csvread( [irfDir, sName, '_IRF_PhatDens_Aggsh1_', num2str(pp), '.csv'], 1, 0);
end

densvalues_ss = PhatDens_IRF(1,:); % steady state density

PhatDens_IRF_5th = zeros(H+1,xn);
PhatDens_IRF_50th = zeros(H+1,xn);
PhatDens_IRF_95th = zeros(H+1,xn);

for hh=2:(H+1)
    PhatDens_IRF_5th(hh,:) =quantile(squeeze(PhatDens_IRF_uncertainty(hh,:,:))-repmat(densvalues_ss',1,Ndraws), 0.05, 2); 
    PhatDens_IRF_50th(hh,:)=quantile(squeeze(PhatDens_IRF_uncertainty(hh,:,:))-repmat(densvalues_ss',1,Ndraws), 0.5, 2); 
    PhatDens_IRF_95th(hh,:)=quantile(squeeze(PhatDens_IRF_uncertainty(hh,:,:))-repmat(densvalues_ss',1,Ndraws), 0.95, 2); 
end

horizons = [5 15 25 50];

%--------------------------------------------------------------------------
% Density Differentials 
%--------------------------------------------------------------------------

for pp = 1:length(horizons)

    hh = horizons(pp);
    
    figure(pp);clf;
    set(figure(pp),'PaperType','usletter','PaperOrientation','Landscape','PaperPosition',[0.1 0.1 11 8.5]);
    plot(xgrid, zeros(size(xgrid)), 'Color', 'k', 'LineStyle', '-', 'LineWidth',2)
    hold on
    plot(xgrid, PhatDens_IRF_5th(hh+1,:), 'Color', 'b', 'LineStyle', ':', 'LineWidth',2)
    hold on
    plot(xgrid, PhatDens_IRF_50th(hh+1,:), 'Color', 'b', 'LineStyle', '-', 'LineWidth',4)
    hold on
    plot(xgrid, PhatDens_IRF_95th(hh+1,:), 'Color', 'b', 'LineStyle', ':', 'LineWidth',2)
    hold on
    plot(pop_Dens_IRF(:,1), pop_Dens_IRF(:,hh+1+1)-pop_Dens_IRF(:,1+1), 'Color', 'k', 'LineStyle', '-', 'LineWidth',4)
    set(gca,'FontSize',50)
    axis([0 2.6 -0.05 0.05])

    sNameDir = ['fVAR', nfVARSpec, '_', nMod, nModSpec, '_MCMC', nMCMCSpec];
    figsaveDir = [pwd, '/', 'Figures/Para', nDataSel ,'/', sNameDir,'/'];
    [~, ~, ~] = mkdir(figsaveDir);

    sNameFile = ['DensityDiff_hh',num2str(hh),'.pdf'];    

    saveas(figure(pp), [figsaveDir sNameFile] );
    %close all
    
end

